model{
  for (i in 1:n) {
    for (j in 1:m) {
      for (k in 1:t) {
        y[i, j, k] ~ dnorm(mu[i, j, k], inv.phi2[j])
        mu[i, j, k] <- beta.raw[j, k] * theta.raw[i, k]
      }
    }
  }
  for (i in 1:n) {
    theta.raw[i, 1] ~ dnorm(0, inv.xi2)
    for (k in 2:t) {
      theta.raw[i, k] ~ dnorm(theta.raw[i, k - 1], inv.tau2.raw)
    }
    for (k in 1:t) {
      theta[i, k] <- (step(beta.raw[1, k]) - (1 - step(beta.raw[1, k]))) *
        (theta.raw[i, k] * inv.xi)
    }
  }
  for (j in 1:m) {
    beta.raw[j, 1] ~ dnorm(0, 1)
    for (k in 2:t) {
      beta.raw[j, k] ~ dnorm(beta.raw[j, k - 1], inv.omega2.raw)
    }
    for (k in 1:t) {
      beta[j, k] <- (step(beta.raw[1, k]) - (1 - step(beta.raw[1, k]))) *
        (beta.raw[j, k] * xi)
    }
    inv.phi2[j] ~ dgamma(1, 0.2)
  }
  phi <- 1 / sqrt(inv.phi2)
  inv.xi2 ~ dgamma(0.5, 0.5)
  inv.xi <- sqrt(inv.xi2)
  xi <- 1 / inv.xi
  kappa ~ dunif(0, 1)
  lambda ~ dunif(0, 1)
  inv.tau2.raw <- inv.xi2 / (kappa * kappa)
  inv.omega2.raw <- 1 / (lambda * lambda)
  tau <- kappa
  omega <- lambda * xi
}